#ifndef DW_MATHS_MATRIX_INL
#	define DW_MATHS_MATRIX_INL

namespace dw {
	namespace maths {

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >::Matrix()
{
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >::Matrix(const T* m)
{
	dwAssert(m != NULL, "Invalid parameter: m is null");
	memcpy(this->m, m, sizeof (T) * M * N);
}

template< typename T, unsigned M, unsigned N >
inline const T& Matrix< T, M, N >::operator[](const unsigned i) const
{ 
	dwAssert(i < M * N, "Invalid parameter: i is out of range");
	return this->m[i];
}

template< typename T, unsigned M, unsigned N >
inline const T& Matrix< T, M, N >::operator()(const unsigned i, const unsigned j) const 
{ 
	dwAssert((j < M) && (i < N), "Invalid parameter: i or j is out of range");
	return this->m[j * N + i];
}

template< typename T, unsigned M, unsigned N >
inline T& Matrix< T, M, N >::operator[](const unsigned i) 
{ 
	dwAssert(i < M * N, "Invalid parameter: i is out of range");
	return this->m[i]; 
}

template< typename T, unsigned M, unsigned N >
inline T& Matrix< T, M, N >::operator()(const unsigned i, const unsigned j)
{ 
	dwAssert((j < M) && (i < N), "Invalid parameter: i or j is out of range");
	return this->m[j * N + i];
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >& Matrix< T, M, N >::operator=(const Matrix< T, M, N >& rhs)
{
	memcpy(&this->m, &rhs.m, sizeof (T) * M * N);
	return *this;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >& Matrix< T, M, N >::operator+=(const Matrix< T, M, N >& rhs)
{	
	for (unsigned i = 0; i < M * N; ++i)
		this->m[i] += rhs.m[i];
	return *this;											
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >& Matrix< T, M, N >::operator-=(const Matrix< T, M, N >& rhs)
{
	for (unsigned i = 0; i < M * N; ++i)
		this->m[i] -= rhs.m[i];
	return *this;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >& Matrix< T, M, N >::operator*=(const T& rhs)
{
	for (unsigned i = 0; i < M * N; ++i)
		this->m[i] *= rhs;
	return *this;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >& Matrix< T, M, N >::operator/=(const T& rhs)
{
	for (unsigned i = 0; i < M * N; ++i)
		this->m[i] /= rhs;
	return *this;
}

template< typename T, unsigned M, unsigned N >
inline bool operator==(const Matrix< T, M, N >& lhs, const Matrix< T, M, N >& rhs)
{
	for (unsigned i = 0; i < M * N; ++i)
		if (lhs.m[i] != rhs.m[i])
			return false;
	return true;
}

template< typename T, unsigned M, unsigned N >
inline bool operator!=(const Matrix< T, M, N >& lhs, const Matrix< T, M, N >& rhs)
{
	return !(lhs == rhs);
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N >& Matrix< T, M, N >::reset()
{
	memset(this->m, 0, sizeof (T) * M * N);
	return *this;
}

template< typename T, unsigned M, unsigned N >
const Vector< T, N >&	Matrix< T, M, N >::getRow(const unsigned i) const
{ 
	dwAssert(i < M, "Invalid parameter: i is out of range");
	return this->v[i];
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N > operator+(const Matrix< T, M, N >& lhs, const Matrix< T, M, N >& rhs)
{
	Matrix< T, M, N > nrv(lhs);
	nrv += rhs;
	return nrv;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N > operator-(const Matrix< T, M, N >& rhs)
{
	return Matrix< T, M, N >::zero -= rhs;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N > operator-(const Matrix< T, M, N >& lhs, const Matrix< T, M, N >& rhs)
{
	Matrix< T, M, N > nrv(lhs);
	nrv -= rhs;
	return nrv;
}

template< typename T, unsigned M, unsigned N, unsigned P >
inline Matrix< T, M, P > operator*(const Matrix< T, M, N >& lhs, const Matrix< T, N, P >& rhs)
{
	Matrix< T, M, P > nrv;

	for (unsigned j = 0; j < M; ++j)
		for (unsigned i = 0; i < P; ++i)
			for (unsigned k = 0; k < N; ++k)
				nrv(i, j) += lhs(i, k) * rhs(k, j);
	return nrv;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N > operator* (const Matrix< T, M, N > lhs, const T& rhs)
{
	Matrix< T, M, N > nrv(lhs);
	nrv *= rhs;
	return nrv;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M, N > operator* (const T& lhs, const Matrix< T, M, N > rhs)
{ 
	Matrix< T, M, N > nrv(rhs);
	nrv *= lhs;
	return nrv;
}

template< typename T, unsigned M, unsigned N >
inline Matrix< T, M - 1, N - 1 > submat(const Matrix< T, M, N >& m, const unsigned col, const unsigned row)
{
	unsigned l, k;
	Matrix< T, M - 1, N - 1 > nrv;

	l = 0;
	for (unsigned j = 0; j < M; ++j)
	{
		k = 0;
		if (j != row)
		{
			for (unsigned i = 0; i < N; ++i)
			if (i != col)
			{
				nrv(k, l) = m(i, j);
				k++;
			}
			l++;
		}
	}
	return nrv;
}

template< typename T, unsigned M, unsigned N >
inline ::std::ostream& operator<<(::std::ostream& os, const Matrix< T, M, N >& m)
{
	for (unsigned j = 0; j < M; ++j)
	{
		os << "| ";
		for (unsigned i = 0; i < N ; ++i)
		{	
			os << ::std::setw(8) << ::std::setprecision(3) << m(i, j);
			if (i < N - 1)
				os << ", ";
		}
		os << " |";
		if (j < M - 1)
			os << "\n";
	}
	return os;
}



// Matrix< T, N, N > partial specialization /////////////////////////////////



template< typename T, unsigned N >
inline Matrix< T, N, N >::Matrix()
{
}

template< typename T, unsigned N >
inline Matrix< T, N, N >::Matrix(const T* m)
{
	dwAssert(m != NULL, "Invalid parameter: m is NULL");
	memcpy(this->m, m, sizeof (T) * N * N);
}

template< typename T, unsigned N >
inline const T& Matrix< T, N, N >::operator[](const unsigned i) const
{ 
	dwAssert(i < N * N, "Invalid parameter: i is out of range");
	return this->m[i];
}

template< typename T, unsigned N >
inline const T& Matrix< T, N, N >::operator()(const unsigned i, const unsigned j) const 
{ 
	dwAssert((j < M) && (i < N), "Invalid parameter: i or j is out of range");
	return this->m[j * N + i];
}

template< typename T, unsigned N >
inline T& Matrix< T, N, N >::operator[](const unsigned i) 
{ 
	dwAssert(i < N * N, "Invalid parameter: i is out of range");
	return this->m[i]; 
}

template< typename T, unsigned N >
inline T& Matrix< T, N, N >::operator()(const unsigned i, const unsigned j)
{ 
	dwAssert((j < N) && (i < N), "Invalid parameter: i or j is out of range");
	return this->m[j * N + i];
}

template< typename T, unsigned N >
inline Matrix< T, N, N >& Matrix< T, N, N >::operator=(const Matrix< T, N, N >& rhs)
{
	memcpy(&this->m, &rhs.m, sizeof (T) * N * N);
	return *this;
}

template< typename T, unsigned N >
inline Matrix< T, N, N >& Matrix< T, N, N >::operator+=(const Matrix< T, N, N >& rhs)
{	
	for (unsigned i = 0; i < N * N; ++i)
		this->m[i] += rhs.m[i];
	return *this;											
}

template< typename T, unsigned N >
inline Matrix< T, N, N >& Matrix< T, N, N >::operator-=(const Matrix< T, N, N >& rhs)
{
	for (unsigned i = 0; i < N * N; ++i)
		this->m[i] -= rhs.m[i];
	return *this;
}

template< typename T, unsigned N >
inline Matrix< T, N, N >& Matrix< T, N, N >::operator*=(const T& rhs)
{
	for (unsigned i = 0; i < M * N; ++i)
		this->m[i] *= rhs;
	return *this;
}

template< typename T, unsigned N >
inline Matrix< T, N, N >& Matrix< T, N, N >::operator/=(const T& rhs)
{
	for (unsigned i = 0; i < N * N; ++i)
		this->m[i] /= rhs;
	return *this;
}

template< typename T, unsigned N >
inline Matrix< T, N, N >& Matrix< T, N, N >::reset()
{
	memset(this->m, 0, sizeof (T) * N * N);
	return *this;
}

template< typename T, unsigned N >
inline const Vector< T, N >&	Matrix< T, N, N >::getRow(const unsigned i) const
{ 
	dwAssert(i < N, "Invalid parameter: i is out of range");
	return this->v[i];
}

template< typename T, unsigned N >
inline T Matrix< T, N, N >::det(const Matrix< T, N, N >& m)
{
	T nrv(0);

	for (unsigned i = 0; i < N; ++i)
		nrv += static_cast< T >(pow(static_cast< T >(-1), static_cast< float >(i))) * m(i, 0) * Matrix< T, N - 1, N - 1 >::det(submat(m, i, 0));
	return nrv;
}

template< typename T, unsigned N >
inline T Matrix< T, N, N >::trace(const Matrix< T, N, N >& m)
{
	T nrv(0);

	for (unsigned i = 0; i < N; ++i)
		nrv += m(i, i);
	return nrv;
}

template< typename T, unsigned N >
inline Matrix< T, N, N > Matrix< T, N, N >::transpose(const Matrix< T, N, N >& m)
{
	Matrix< T, N, N > nrv;
	
	for (unsigned j = 0; j < N; ++j)
		for (unsigned i = 0; i < N; ++i)
			nrv(i, j) = m(j, i);
	return nrv;
}

template< typename T, unsigned N >
inline Matrix< T, N, N > Matrix< T, N, N >::cofactor(const Matrix< T, N, N >& m)
{
	Matrix< T, N, N > nrv;

	for (unsigned j = 0; j < N; ++j)
		for (unsigned i = 0; i < N; ++i)
			nrv(i, j) = std::pow(-1, i + j) * Matrix< T, N - 1, N - 1 >::det(submat(m, i, j));
	return nrv;
}

template< typename T, unsigned N >
inline Matrix< T, N, N > Matrix< T, N, N >::adjoint(const Matrix< T, N, N >& m)
{
	return transpose(cofactor(m));
}

template< typename T, unsigned N >
inline Matrix< T, N, N > Matrix< T, N, N >::inverse(const Matrix< T, N, N >& m)
{ 
	return 1 / det(m) * adjoint(m);
}



// Matrix< T, 2, 2 > partial specialization /////////////////////////////////



template< typename T > const Matrix< T, 2, 2 > Matrix< T, 2, 2 >::zero(0, 0, 0, 0);

template< typename T > const Matrix< T, 2, 2 > Matrix< T, 2, 2 >::identity(1, 0, 0, 1);

template< typename T >
inline Matrix< T, 2, 2 >::Matrix()
{
}

template< typename T >
inline Matrix< T, 2, 2 >::Matrix(const T* m)
{
	dwAssert(m != NULL, "Invalid parameter: m is NULL");
	memcpy(this->m, m, sizeof (T) * 4);
}

template< typename T >
inline Matrix< T, 2, 2 >::Matrix(
	const T& m00, const T& m01, 
	const T& m10, const T& m11) :
	m00(m00), m01(m01),
	m10(m10), m11(m11)
{
}

template< typename T >
inline Matrix< T, 2, 2 >::Matrix(const Vector< T, 2 >& v0, const Vector< T, 2 >& v1) :
	v0(v0), v1(v1)
{
}

template< typename T >
inline const T& Matrix< T, 2, 2 >::operator[](const unsigned i) const
{ 
	dwAssert(i < 4, "Invalid parameter: i is out of range");
	return this->m[i];
}

template< typename T >
inline const T& Matrix< T, 2, 2 >::operator()(const unsigned i, const unsigned j) const 
{ 
	dwAssert((j < 2) && (i < 2), "Invalid parameter: i or j is out of range");
	return this->m[j * 2 + i];
}

template< typename T >
inline T& Matrix< T, 2, 2 >::operator[](const unsigned i) 
{ 
	dwAssert(i < 4, "Invalid parameter: i is out of range");
	return this->m[i]; 
}

template< typename T >
inline T& Matrix< T, 2, 2 >::operator()(const unsigned i, const unsigned j)
{ 
	dwAssert((j < 2) && (i < 2), "Invalid parameter: i or j is out of range");
	return this->m[j * 2 + i];
}

template< typename T >
inline Matrix< T, 2, 2 >& Matrix< T, 2, 2 >::operator=(const Matrix< T, 2, 2 >& rhs)
{
	memcpy(&this->m, &rhs.m, sizeof (T) * 4);
	return *this;
}

template< typename T >
inline Matrix< T, 2, 2 >& Matrix< T, 2, 2 >::operator+=(const Matrix< T, 2, 2 >& rhs)
{	
	for (unsigned i = 0; i < 4; ++i)
		this->m[i] += rhs.m[i];
	return *this;											
}

template< typename T >
inline Matrix< T, 2, 2 >& Matrix< T, 2, 2 >::operator-=(const Matrix< T, 2, 2 >& rhs)
{
	for (unsigned i = 0; i < 4; ++i)
		this->m[i] -= rhs.m[i];
	return *this;
}

template< typename T >
inline Matrix< T, 2, 2 >& Matrix< T, 2, 2 >::operator*=(const T& rhs)
{
	for (unsigned i = 0; i < 4; ++i)
		this->m[i] *= rhs;
	return *this;
}

template< typename T >
inline Matrix< T, 2, 2 >& Matrix< T, 2, 2 >::operator/=(const T& rhs)
{
	for (unsigned i = 0; i < 2 * 2; ++i)
		this->m[i] /= rhs;
	return *this;
}

template< typename T >
inline const Vector< T, 2 >&	Matrix< T, 2, 2 >::getRow(const unsigned i) const
{ 
	dwAssert(i < 2, "Invalid parameter: i is out of range");
	return this->v[i];
}

template< typename T >
inline T Matrix< T, 2, 2 >::det(const Matrix< T, 2, 2 >& m)
{
	return m[0] * m[3] - m[2] * m[1];
}

template< typename T >
inline T Matrix< T, 2, 2 >::trace(const Matrix< T, 2, 2 >& m)
{
	return m.m00 + m.m11;
}

template< typename T >
inline Matrix< T, 2, 2 > Matrix< T, 2, 2 >::transpose(const Matrix< T, 2, 2 >& m)
{
	return Matrix< T, 2, 2 >(
		m.m00, m.m10,
		m.m01, m.m11);
}

template< typename T >
inline Matrix< T, 2, 2 > Matrix< T, 2, 2 >::cofactor(const Matrix< T, 2, 2 >& m)
{
	Matrix< T, 2, 2 > nrv;

	for (unsigned j = 0; j < 2; ++j)
		for (unsigned i = 0; i < 2; ++i)
			nrv(i, j) = std::pow(static_cast< T >(-1), i + j) * Matrix< T, 1, 1 >::det(submat(m, i, j));
	return nrv;
}

template< typename T >
inline Matrix< T, 2, 2 > Matrix< T, 2, 2 >::adjoint(const Matrix< T, 2, 2 >& m)
{
	return transpose(cofactor(m));
}

template< typename T >
inline Matrix< T, 2, 2 > Matrix< T, 2, 2 >::inverse(const Matrix< T, 2, 2 >& m)
{ 
	return 1 / det(m) * adjoint(m);
}



// Matrix< T, 3, 3 > partial specialization /////////////////////////////////



template< typename T > const Matrix< T, 3, 3 > Matrix< T, 3, 3 >::zero(0, 0, 0, 0, 0, 0, 0, 0, 0);

template< typename T > const Matrix< T, 3, 3 > Matrix< T, 3, 3 >::identity(1, 0, 0, 0, 1, 0, 0, 0, 1);

template< typename T >
inline Matrix< T, 3, 3 >::Matrix()
{
}

template< typename T >
inline Matrix< T, 3, 3 >::Matrix(const T* m)
{
	dwAssert(m != NULL, "Invalid parameter: m is NULL");
	memcpy(this->m, m, sizeof (T) * 9);
}

template< typename T >
inline Matrix< T, 3, 3 >::Matrix(
	const T& m00, const T& m01, const T& m02, 
	const T& m10, const T& m11, const T& m12, 
	const T& m20, const T& m21, const T& m22) :
	right(m00, m01, m02),
	up(m10, m11, m12),
	at(m20, m21, m22)
{
}

template< typename T >
inline Matrix< T, 3, 3 >::Matrix(const Vector< T, 3 >& v0, const Vector< T, 3 >& v1, const Vector< T, 3 >& v2) :
	right(v0), up(v1), at(v2)
{
}

template< typename T >
inline const T& Matrix< T, 3, 3 >::operator[](const unsigned i) const
{ 
	dwAssert(i < 9, "Invalid parameter: i is out of range");
	return this->m[i];
}

template< typename T >
inline const T& Matrix< T, 3, 3 >::operator()(const unsigned i, const unsigned j) const 
{ 
	dwAssert((j < 3) && (i < 3), "Invalid parameter: i or j is out of range");
	return this->m[j * 3 + i];
}

template< typename T >
inline T& Matrix< T, 3, 3 >::operator[](const unsigned i) 
{ 
	dwAssert(i < 9, "Invalid parameter: i is out of range");
	return this->m[i]; 
}

template< typename T >
inline T& Matrix< T, 3, 3 >::operator()(const unsigned i, const unsigned j)
{ 
	dwAssert((j < 3) && (i < 3), "Invalid parameter: i or j is out of range");
	return this->m[j * 3 + i];
}

template< typename T >
inline Matrix< T, 3, 3 >& Matrix< T, 3, 3 >::operator=(const Matrix< T, 3, 3 >& rhs)
{
	memcpy(&this->m, &rhs.m, sizeof (T) * 9);
	return *this;
}

template< typename T >
inline Matrix< T, 3, 3 >& Matrix< T, 3, 3 >::operator+=(const Matrix< T, 3, 3 >& rhs)
{	
	for (unsigned i = 0; i < 9; ++i)
		this->m[i] += rhs.m[i];
	return *this;											
}

template< typename T >
inline Matrix< T, 3, 3 >& Matrix< T, 3, 3 >::operator-=(const Matrix< T, 3, 3 >& rhs)
{
	for (unsigned i = 0; i < 9; ++i)
		this->m[i] -= rhs.m[i];
	return *this;
}

template< typename T >
inline Matrix< T, 3, 3 >& Matrix< T, 3, 3 >::operator*=(const T& rhs)
{
	for (unsigned i = 0; i < 9; ++i)
		this->m[i] *= rhs;
	return *this;
}

template< typename T >
inline Matrix< T, 3, 3 >& Matrix< T, 3, 3 >::operator/=(const T& rhs)
{
	for (unsigned i = 0; i < 9; ++i)
		this->m[i] /= rhs;
	return *this;
}

template< typename T >
inline const Vector< T, 3 >&	Matrix< T, 3, 3 >::getRow(const unsigned i) const
{ 
	assert(i < 3);
	return this->v[i];
}

template< typename T >
inline T Matrix< T, 3, 3 >::det(const Matrix< T, 3, 3 >& m)
{
	return 
		m[0] * m[4] * m[8] - m[0] * m[7] * m[5] 
	- m[3] * m[1] * m[8] + m[3] * m[7] * m[2] 
	+ m[6] * m[1] * m[5] - m[6] * m[4] * m[2];
}

template< typename T >
inline T Matrix< T, 3, 3 >::trace(const Matrix< T, 3, 3 >& m)
{
	return m.m00 + m.m11 + m.m22;
}

template< typename T >
inline Matrix< T, 3, 3 > Matrix< T, 3, 3 >::transpose(const Matrix< T, 3, 3 >& m)
{
	return Matrix< T, 3, 3 >(
		m.m00, m.m10, m.m20,
		m.m01, m.m11, m.m21,
		m.m02, m.m12, m.m22);
}

template< typename T >
inline Matrix< T, 3, 3 > Matrix< T, 3, 3 >::cofactor(const Matrix< T, 3, 3 >& m)
{
	Matrix< T, 3, 3 > nrv;

	for (dw::uint32 j = 0; j < 3; ++j)
		for (dw::uint32 i = 0; i < 3; ++i)
			nrv(i, j) = std::pow(static_cast< T >(-1), i + j) * Matrix< T, 2, 2 >::det(submat(m, i, j));
	return nrv;
}

template< typename T >
inline Matrix< T, 3, 3 > Matrix< T, 3, 3 >::adjoint(const Matrix< T, 3, 3 >& m)
{
	return transpose(cofactor(m));
}

template< typename T >
inline Matrix< T, 3, 3 > Matrix< T, 3, 3 >::inverse(const Matrix< T, 3, 3 >& m)
{ 
	return 1 / det(m) * adjoint(m);
}

template< typename T >
inline Matrix< T, 3, 3 >	Matrix< T, 3, 3 >::rotation(const Vector< T, 3 >& axis, const T& angle)
{
	Quaternion< T > q(axis, angle);

	return q.toMatrix();
}

template< typename T >
inline Matrix< T, 3, 3 >	Matrix< T, 3, 3 >::scaling(const Vector< T, 3 >& v)
{
	Matrix< T, 3, 3 > nrv(Matrix< T, 3, 3 >::zero);

	nrv.right.x = v[0];
	nrv.up.y = v[1];
	nrv.at.z = v[2];
	return nrv;
}

template< typename T >
inline Matrix< T, 3, 3 > operator*(const Matrix< T, 3, 3 >& lhs, const Matrix< T, 3, 3 >& rhs)
{
	Matrix< T, 3, 3 > nrv;

	nrv[0] = lhs[0] * rhs[0] + lhs[1] * rhs[3] + lhs[2] * rhs[6];
	nrv[1] = lhs[0] * rhs[1] + lhs[1] * rhs[4] + lhs[2] * rhs[7];
	nrv[2] = lhs[0] * rhs[2] + lhs[1] * rhs[5] + lhs[2] * rhs[8];
	nrv[3] = lhs[3] * rhs[0] + lhs[4] * rhs[3] + lhs[5] * rhs[6];
	nrv[4] = lhs[3] * rhs[1] + lhs[4] * rhs[4] + lhs[5] * rhs[7];
	nrv[5] = lhs[3] * rhs[2] + lhs[4] * rhs[5] + lhs[5] * rhs[8];
	nrv[6] = lhs[6] * rhs[0] + lhs[7] * rhs[3] + lhs[8] * rhs[6];
	nrv[7] = lhs[6] * rhs[1] + lhs[7] * rhs[4] + lhs[8] * rhs[7];
	nrv[8] = lhs[6] * rhs[2] + lhs[7] * rhs[5] + lhs[8] * rhs[8];
	return nrv;
}

//==============================================================================
// Matrix<4, 4, T>
//==============================================================================

//------------------------------------------------------------------------------
// Constants
//------------------------------------------------------------------------------

template< typename T >
/*static*/const Matrix< T, 4, 4 > Matrix< T, 4, 4 >::zero(0, 0, 0, 0,
														  0, 0, 0, 0,
														  0, 0, 0, 0,
														  0, 0, 0, 0);

template< typename T >
/*static*/ const Matrix< T, 4, 4 > Matrix< T, 4, 4 >::identity(1, 0, 0, 0,
															   0, 1, 0, 0,
															   0, 0, 1, 0,
															   0, 0, 0, 1);

//------------------------------------------------------------------------------
// Construction & Destruction
//------------------------------------------------------------------------------
/*
//------------------------------------------------------------------------------
template< typename T >
inline Matrix< T, 4, 4 >::Matrix()
:	right(Vector< T, 4 >::unitX)
,	up(Vector< T, 4 >::unitY)
,	at(Vector< T, 4 >::unitZ)
,	pos(Vector< T, 4 >::unitW)
{
}
*/
//------------------------------------------------------------------------------
template< typename T >
inline Matrix< T, 4, 4 >::Matrix(const Vector< T, 4 >& pRight, const Vector< T, 4 >& pUp, const Vector< T, 4 >& pAt, const Vector< T, 4 >& pPos)
:	right(pRight)
,	up(pUp)
,	at(pAt)
,	pos(pPos)
{
}

//------------------------------------------------------------------------------
template< typename T >
inline Matrix< T, 4, 4 >::Matrix(const Matrix< T, 4, 4 >& ref /*= Identity*/)
:	right(ref.right)
,	up(ref.up)
,	at(ref.at)
,	pos(ref.pos)
{
}

//------------------------------------------------------------------------------
template< typename T >
inline Matrix< T, 4, 4 >::Matrix(const T* m)
:	right(&m[0])
,	up(&m[4])
,	at(&m[8])
,	pos(&m[12])
{
	dwAssert(m != NULL, "Invalid parameter: m is NULL");
}

//------------------------------------------------------------------------------
template< typename T >
inline Matrix< T, 4, 4 >::Matrix(const T& m00, const T& m01, const T& m02, const T& m03,
							   const T& m10, const T& m11, const T& m12, const T& m13,
							   const T& m20, const T& m21, const T& m22, const T& m23,
							   const T& m30, const T& m31, const T& m32, const T& m33)
:	right(m00, m01, m02, m03)
,	up(m10, m11, m12, m13)
,	at(m20, m21, m22, m23)
,	pos(m30, m31, m32, m33)
{
}

//------------------------------------------------------------------------------
template< typename T >
inline Matrix< T, 4, 4 >::Matrix(const Matrix< T, 3, 3 >& rot, const Vector< T, 3 >& pPos)
:	right(rot.right, 0)
,	up(rot.up, 0)
,	at(rot.at, 0)
,	pos(pPos, 1)
{
}

//------------------------------------------------------------------------------
template< typename T >
inline const T& Matrix< T, 4, 4 >::operator[](const unsigned i) const
{ 
	dwAssert(i < 16, "Invalid parameter: i is out of range");
	return this->m[i];
}

template< typename T >
inline const T& Matrix< T, 4, 4 >::operator()(const unsigned i, const unsigned j) const 
{ 
	dwAssert((j < 4) && (i < 4), "Invalid parameter: i or j is out of range");
	return this->m[j * 4 + i];
}

template< typename T >
inline T& Matrix< T, 4, 4 >::operator[](const unsigned i) 
{ 
	dwAssert(i < 16, "Invalid parameter: i is out of range");
	return this->m[i]; 
}

template< typename T >
inline T& Matrix< T, 4, 4 >::operator()(const unsigned i, const unsigned j)
{ 
	dwAssert((j < 4) && (i < 4), "Invalid parameter: i or j is out of range");
	return this->m[j * 4 + i];
}

template< typename T >
inline Matrix< T, 4, 4 >& Matrix< T, 4, 4 >::operator=(const Matrix< T, 4, 4 >& rhs)
{
	memcpy(&this->m, &rhs.m, sizeof (T) * 4 * 4);
	return *this;
}

template< typename T >
inline Matrix< T, 4, 4 >& Matrix< T, 4, 4 >::operator+=(const Matrix< T, 4, 4 >& rhs)
{	
	for (unsigned i = 0; i < 16; ++i)
		this->m[i] += rhs.m[i];
	return *this;											
}

template< typename T >
inline Matrix< T, 4, 4 >& Matrix< T, 4, 4 >::operator-=(const Matrix< T, 4, 4 >& rhs)
{
	for (unsigned i = 0; i < 16; ++i)
		this->m[i] -= rhs.m[i];
	return *this;
}

template< typename T >
inline Matrix< T, 4, 4 >& Matrix< T, 4, 4 >::operator*=(const T& rhs)
{
	for (unsigned i = 0; i < 16; ++i)
		this->m[i] *= rhs;
	return *this;
}

template< typename T >
inline Matrix< T, 4, 4 >& Matrix< T, 4, 4 >::operator/=(const T& rhs)
{
	for (unsigned i = 0; i < 16; ++i)
		this->m[i] /= rhs;
	return *this;
}

template< typename T >
inline const Vector< T, 4 >&	Matrix< T, 4, 4 >::getRow(const unsigned i) const
{ 
	dwAssert(i < 4, "Invalid parameter: i is out of range");
	return this->v[i];
}

template< typename T >
/*static*/ inline T Matrix< T, 4, 4 >::det(const Matrix< T, 4, 4 >& m)
{
	const T a = m.at.z * m.pos.w - m.at.w * m.pos.z; //
	const T b = m.at.y * m.pos.w - m.at.w * m.pos.y;
	const T c = m.at.y * m.pos.z - m.at.z * m.pos.y; //
	const T d = m.at.x * m.pos.w - m.at.w * m.pos.x;
	const T e = m.at.x * m.pos.z - m.at.z * m.pos.x;
	const T f = m.at.x * m.pos.y - m.at.y * m.pos.x; //
	const T g = m.at.x * m.pos.y - m.at.y * m.pos.w;
	
	return m.right.x * (m.up.y * a - m.up.z * b + m.up.w * c) -
		   m.right.y * (m.up.x * a - m.up.z * d + m.up.w * e) +
           m.right.z * (m.up.x * b - m.up.y * d + m.up.w * f) -
           m.right.w * (m.up.x * c - m.up.y * e + m.up.z * g);
}

template< typename T >
inline Matrix< T, 4, 4 > Matrix< T, 4, 4 >::transpose(const Matrix< T, 4, 4 >& m)
{
	return Matrix< T, 4, 4 >(
		m[0, 0], m[1, 0], m[2, 0], m[3, 0],
		m[0, 1], m[1, 1], m[2, 1], m[3, 1],
		m[0, 2], m[1, 2], m[2, 2], m[3, 2],
		m[0, 3], m[1, 3], m[2, 3], m[3, 3]);
}

template< typename T >
inline Matrix< T, 4, 4 > Matrix< T, 4, 4 >::inverse(const Matrix< T, 4, 4 >& m)
{
	const T determinant = det(m);

	dwAssert(determinant != 0, "");

	const T invDet = 1 / determinant;

	Matrix< T, 4, 4 > invMtx;

	invMtx.right.x =  (m.up.y    * (m.at.z * m.pos.w - m.pos.z * m.at.w) - m.up.z * (m.at.y * m.pos.w - m.at.w * m.pos.y) + m.up.w * (m.at.y * m.pos.z - m.at.z *  m.pos.y)) * invDet;
    invMtx.right.y = -(m.right.y * (m.at.z * m.pos.w - m.pos.z * m.at.w) - m.at.y * (m.right.z * m.pos.w - m.pos.z * m.right.w) + m.pos.y * (m.right.z * m.at.w - m.at.z *  m.right.w)) * invDet;
    invMtx.right.z =  (m.right.y * (m.up.z * m.pos.w - m.pos.z * m.up.w) - m.up.y * (m.right.z * m.pos.w - m.pos.z * m.right.w) + m.pos.y * (m.right.z * m.up.w - m.up.z *  m.right.w)) * invDet;
    invMtx.right.w = -(m.right.y * (m.up.z * m.at.w - m.at.z * m.up.w) - m.up.y * (m.right.z * m.at.w - m.at.z * m.right.w) + m.at.y * (m.right.z * m.up.w - m.up.z *  m.right.w)) * invDet;
        
	invMtx.up.x = -(m.up.x    * (m.at.z * m.pos.w - m.at.w * m.pos.z) - m.up.z * (m.at.x * m.pos.w - m.at.w * m.pos.x) + m.up.w * (m.at.x * m.pos.z - m.at.z *  m.pos.x)) * invDet;
    invMtx.up.y =  (m.right.x * (m.at.z * m.pos.w - m.pos.z * m.at.w) - m.at.x * (m.right.z * m.pos.w - m.pos.z * m.right.w) + m.pos.x * (m.right.z * m.at.w - m.at.z *  m.right.w)) * invDet;
    invMtx.up.z = -(m.right.x * (m.up.z * m.pos.w - m.pos.z * m.up.w) - m.up.x * (m.right.z * m.pos.w - m.pos.z * m.right.w) + m.pos.x * (m.right.z * m.up.w - m.up.z *  m.right.w)) * invDet;
    invMtx.up.w =  (m.right.x * (m.up.z * m.at.w - m.at.z * m.up.w) - m.up.x * (m.right.z * m.at.w - m.at.z * m.right.w) + m.at.x * (m.right.z * m.up.w - m.up.z *  m.right.w)) * invDet;
        
	invMtx.at.x =  (m.up.x    * (m.at.y * m.pos.w - m.at.w * m.pos.y) - m.up.y * (m.at.x * m.pos.w - m.at.w * m.pos.x) + m.up.w * (m.at.x * m.pos.y - m.at.y *  m.pos.x)) * invDet;
    invMtx.at.y = -(m.right.x * (m.at.y * m.pos.w - m.pos.y * m.at.w) - m.at.x * (m.right.y * m.pos.w - m.pos.y * m.right.w) + m.pos.x * (m.right.y * m.at.w - m.at.y *  m.right.w)) * invDet;
    invMtx.at.z =  (m.right.x * (m.up.y * m.pos.w - m.pos.y * m.up.w) - m.up.x * (m.right.y * m.pos.w - m.pos.y * m.right.w) + m.pos.x * (m.right.y * m.up.w - m.up.y *  m.right.w)) * invDet;
    invMtx.at.w = -(m.right.x * (m.up.y * m.at.w - m.at.y * m.up.w) - m.up.x * (m.right.y * m.at.w - m.at.y * m.right.w) + m.at.x * (m.right.y * m.up.w - m.up.y *  m.right.w)) * invDet;
        
	invMtx.pos.x = -(m.up.x    * (m.at.y * m.pos.z - m.at.z * m.pos.y) - m.up.y * (m.at.x * m.pos.z - m.at.z * m.pos.x) + m.up.z * (m.at.x * m.pos.y - m.at.y *  m.pos.x)) * invDet;
    invMtx.pos.y =  (m.right.x * (m.at.y * m.pos.z - m.pos.y * m.at.z) - m.at.x * (m.right.y * m.pos.z - m.pos.y * m.right.z) + m.pos.x * (m.right.y * m.at.z - m.at.y *  m.right.z)) * invDet;
    invMtx.pos.z = -(m.right.x * (m.up.y * m.pos.z - m.pos.y * m.up.z) - m.up.x * (m.right.y * m.pos.z - m.pos.y * m.right.z) + m.pos.x * (m.right.y * m.up.z - m.up.y *  m.right.z)) * invDet;
	invMtx.pos.w =  (m.right.x * (m.up.y * m.at.z - m.at.y * m.up.z) - m.up.x * (m.right.y * m.at.z - m.at.y * m.right.z) + m.at.x * (m.right.y * m.up.z - m.up.y *  m.right.z)) * invDet;

	return invMtx;
}

template< typename T >
inline Matrix< T, 4, 4 >	Matrix< T, 4, 4 >::translation(const Vector< T, 3 >& v)
{
	Matrix< T, 4, 4 > nrv(Matrix< T, 4, 4 >::identity);

	nrv.pos.x = v[0];
	nrv.pos.y = v[1];
	nrv.pos.z = v[2];

	return nrv;
}

template< typename T >
/*static*/ Matrix< T, 4, 4 > Matrix< T, 4, 4 >::perspective(T fov, T aspectRatio, T zNear, T zFar)
{
	Matrix< T, 4, 4 > proj;

	float f = 1.0f / ::tanf(fov * 0.5f);
	float coef = 1.0f / (zNear - zFar);

	proj.right.set(f / aspectRatio, 0.0f, 0.0f, 0.0f);
	proj.up.set(0.0f, f, 0.0f, 0.0f);
	proj.at.set(0.0f, 0.0f, (zFar + zNear) * coef, -1.0f);
	proj.pos.set(0.0f, 0.0f, (2.0f * zFar * zNear) * coef, 0.0f);

	return proj;
}

template< typename T >
/*static*/ Matrix< T, 4, 4 > Matrix< T, 4, 4 >::view(const Matrix< T, 4, 4 >& transform)
{
	maths::Matrix4f viewMatrix;

    viewMatrix.m[0] = -transform.right.x;
    viewMatrix.m[1] =  transform.up.x;
    viewMatrix.m[2] = -transform.at.x;
    viewMatrix.m[3] =  0.0f;

    viewMatrix.m[4] = -transform.right.y;
    viewMatrix.m[5] =  transform.up.y;
    viewMatrix.m[6] = -transform.at.y;
    viewMatrix.m[7] =  0.0f;

    viewMatrix.m[8]  = -transform.right.z;
    viewMatrix.m[9]  =  transform.up.z;
    viewMatrix.m[10] = -transform.at.z;
    viewMatrix.m[11] =  0.0f;

	viewMatrix.m[12] =  maths::Vector4f::dot(transform.right, transform.pos);
    viewMatrix.m[13] = -maths::Vector4f::dot(transform.up, transform.pos);
    viewMatrix.m[14] =  maths::Vector4f::dot(transform.at, transform.pos);
    viewMatrix.m[15] =  1.0f;

	return viewMatrix;
}

template< typename T >
inline Matrix< T, 4, 4 > operator*(const Matrix< T, 4, 4 >& lhs, const Matrix< T, 4, 4 >& rhs)
{ 
	Matrix< T, 4, 4 > nrv;

	nrv[0]	= lhs[0]	* rhs[0] + lhs[1]		* rhs[4] + lhs[2]		* rhs[8]	+ lhs[3]	* rhs[12];
	nrv[1]	= lhs[0]	* rhs[1] + lhs[1]		* rhs[5] + lhs[2]		* rhs[9]	+ lhs[3]	* rhs[13];
	nrv[2]	= lhs[0]	* rhs[2] + lhs[1]		* rhs[6] + lhs[2]		* rhs[10] + lhs[3]	* rhs[14];
	nrv[3]	= lhs[0]	* rhs[3] + lhs[1]		* rhs[7] + lhs[2]		* rhs[11] + lhs[3]	* rhs[15];
	nrv[4]	= lhs[4]	* rhs[0] + lhs[5]		* rhs[4] + lhs[6]		* rhs[8]	+ lhs[7]	* rhs[12];
	nrv[5]	= lhs[4]	* rhs[1] + lhs[5]		* rhs[5] + lhs[6]		* rhs[9]	+ lhs[7]	* rhs[13];
	nrv[6]	= lhs[4]	* rhs[2] + lhs[5]		* rhs[6] + lhs[6]		* rhs[10] + lhs[7]	* rhs[14];
	nrv[7]	= lhs[4]	* rhs[3] + lhs[5]		* rhs[7] + lhs[6]		* rhs[11] + lhs[7]	* rhs[15];
	nrv[8]	= lhs[8]	* rhs[0] + lhs[9]		* rhs[4] + lhs[10]	* rhs[8]	+ lhs[11] * rhs[12];
	nrv[9]	= lhs[8]	* rhs[1] + lhs[9]		* rhs[5] + lhs[10]	* rhs[9]	+ lhs[11] * rhs[13];
	nrv[10]	= lhs[8]	* rhs[2] + lhs[9]		* rhs[6] + lhs[10]	* rhs[10] + lhs[11] * rhs[14];
	nrv[11]	= lhs[8]	* rhs[3] + lhs[9]		* rhs[7] + lhs[10]	* rhs[11] + lhs[11] * rhs[15];
	nrv[12]	= lhs[12] * rhs[0] + lhs[13]	* rhs[4] + lhs[14]	* rhs[8]	+ lhs[15] * rhs[12];
	nrv[13]	= lhs[12] * rhs[1] + lhs[13]	* rhs[5] + lhs[14]	* rhs[9]	+ lhs[15] * rhs[13];
	nrv[14]	= lhs[12] * rhs[2] + lhs[13]	* rhs[6] + lhs[14]	* rhs[10]	+ lhs[15] * rhs[14];
	nrv[15]	= lhs[12] * rhs[3] + lhs[13]	* rhs[7] + lhs[14]	* rhs[11] + lhs[15] * rhs[15];
	return nrv;
}

	} // namespace maths
} // namespace dw

#endif // !DW_MATHS_MATRIX_INL
